Numerical simulations on compression behaviors of the laminated shale based on the digital image technology and the discrete element method

As an unconventional reservoir sedimentary rock, shale contains a series of layers and various microstructures that lead to complex mechanical properties, such as the anisotropy of stiffness and strength. This study is directed towards the anisotropy caused by the microstructures of the shale, employing the 2D particle flow code (PFC2D) to explore stiffness, strength, failure mode, and micro-crack evolution. More realistic microstructures and the calibration of microscopic parameters of the shale are reasonably considered through the computed tomography (CT) images and mineral analysis. The corresponding numerical simulation results are fully compared with the experimental results. In what follows, the sensitivity analysis is conducted on the key microscopic parameters and microstructure characteristics in numerical samples with laminated characteristics. The results show that the influence of microscopic parameters of the parallel bonding model on macroscopic parameters is related to the layering angle and the face type, and the microstructures and initial cracks of numerical samples can considerably affect the macroscopic mechanical behaviors of the laminated samples. Next, the effect of confining pressure on the mechanical properties of layered shale is also discussed based on the numerical results. These findings highlight the potential of this approach for applications in micro-scaled models and calibration of microscopic parameters to probe mechanical behaviors of the laminated rock.


Slice images of the microstructures of shale samples
The CT data information of Mancos shale published by Ramos et al. 6,7 in the Digital Rock Portal 33 is used to extract the geometric properties of laminated shale throughout this work.This set of data includes a total of 9 samples with the layering angle between the drilling direction and the bedding direction being 0° (horizontal), 45°, and 90°.All of the samples are cylindrical specimens with the diameter being d = 25 mm and the height being h = 50 mm approximately.The samples are underwent CT scans before and after the triaxial test, thus a total of 18 sets of digital core data can be obtained.
Figure 1 shows the results of 3D reconstruction of the state of each sample before and after the triaxial test.The directions of all sample data were adjusted to make their bedding directions parallel to the y-axis to ensure accurate comparison of pre and post samples.Although this set of data includes 9 samples, most of the sample exist initial defects, which can be observed through the pre sample digital core.These initial defects, which are inevitable due to factors such as rock block sampling or sample cutting, can significantly reduce the stiffness and strength of the sample, and it should be considered when we investigate the impact of sample layering on rock properties.Therefore, three samples different layering angles (red labels in Fig. 1) which are equipped with the least initial defects are selected as the analysis objects.
The slicing planes for the three specimens are along the symmetric axis of the sample and parallel to the x-z plane (as shown in Fig. 2a), to fully reflect the layering characteristic of the sample.By comparing the details of shale minerals in the slice images of pre-and post-test in Fig. 2b-d, it is clearly seen that the angle adjustment of digital cores is accurate and the pre and post sliced images are at the same location.Mancos samples exhibit two major facies in their bedding (Fig. 2a), and both facies contain the same major mineral constituents, but they are in different quantities 7 .The light facies (LF) consists of approximately 52% quartz, 13% clays, 16% calcite, and 11.4% dolomite, whereas the dark facies is comprised of approximately 15% quartz, 46% clays, 5% calcite, and

Image processing of the shale sample slices
The initial sliced image is converted into a grayscale image, which is actually composed of a numerical matrix F M×N and can be expressed as where M and N denote the total number of rows and columns of pixels in the image, respectively, and f(x, y) represents the grayscale value at the pixel point (x, y) in the range [0, 255].
Due to the errors in sample size and the scanning segment length, the size ratio of the initial sliced images of the different samples are inconsistent in Fig. 2b-d.In addition, lower grayscale owing to uneven lighting at the edges of sliced image also lead to errors in mineral recognition based on grayscale.Therefore, initial sliced images should be cut with their center pixel as the center to obtain the analyzed sliced images with the diameter to height ratio being 0.5.Then, moderate gradual adjustment 35 of grayscale at the corresponding edges of each image was applied to eliminate the grayscale errors, and the resulting slice images are shown in Fig. 3.The main features of the images in Fig. 3 are the layering structure between LF and DF.Although there are still some pyrite or feldspar minerals with high grayscale values, as well as micro-cracks or pores with low grayscale values in the image, those components were correspondingly divided into LF or DF due to the low fraction.
(1) Then, the threshold segmentation was applied to divide the regions of LF and DF which can convert the grayscale image into binary images.The Otsu method 36 is the most common and reasonable way to obtain the grayscale value of threshold and as an example, the results of threshold segmentation in the 45° sample slice is shown in Fig. 4a.However, this image still needs further morphological operation 37 in order to be applied to subsequent modeling, and the processes shown in Fig. 4 includes: (1) Remove small spots of LF and DF (by invert the image) to filter out small particles.(2) Closing operation to fill the gaps or defects in the edge of LF or DF whose structures are much smaller than those in the region of layering facies.(3) Opening operation to remove the bulges in the edge of LF or DF.(4) Remove small spots of LF and DF once more.In brief, these operations reduce the complexity of geometry on the premise of maintaining the overall layering properties which can be  www.nature.com/scientificreports/seen in Fig. 4e.By adjusting the corresponding parameters in morphological operation, the area ratio of LF and DF can be controlled to the value of 1:1 approximately.

Establishment of the numerical sample
A 2D discrete element sample with equal proportions of LF and DF was established by using cellular automata simulation 38 .Well-contact disc particles evolved automatically and generated mineral microstructures to simulate the characteristics of rocks based on PFC 2D .For the shale sample in this study, binary sliced images are used to group the generated particles for the subsequent mechanics definition.The numerical sample with the diameter D = 25 mm and the height H = 50 mm is established.It is generally believed that when the particle dimension parameter RES > 10, the influence of particle size on DEM simulation results can be ignored 39 .Based on previous simulation studies on DEM 26,30 , a total of 27,222 discs with particle radius ranging from 0.015 to 0.030 mm were generated which RES > > 500.These parameters were chosen to ensure that the numerical model faithfully represents the heterogeneity and anisotropic properties of the shale samples.Figure 5 shows the generation results of the 45° numerical sample.The overall distribution of mineral facies and the micro-structures in layers are well characterized using the discrete numerical model by comparing Fig. 5a-d, respectively.For convenience of subsequent comparison and discussion, the classical and ideal interbedded model 26 should also be built with the interlayer thickness d = 2.5 mm, which is shown in Fig. 5e.

Linear parallel bond model
The particle DEM could be used to simulate the mineral particles of rock as discs in 2D that are in contact with each other, and it can categorize the properties of discs and contacts to simulate microscopic rock minerals 21 .The linear parallel-bond model (Fig. 6) proposed by Potyondy 19,40 is appropriate for simulating the micromechanical properties of cemented materials, such as rocks.The parallel-bond (PB) model consists of the linear contact element and the parallel-bond element as shown in Fig. 6.The linear contact element (red line in Fig. 6) is equivalent to the linear contact model which can only transmit pressure and friction-sliding and cannot resist tension and rotation.The parallel-bond element (blue line in Fig. 6) is called a parallel bond which acts in parallel with the linear contact element and can resist tension and rotation.When the parallel bond does not work, the linear parallel-bond model will degenerate into the linear contact model.
There are several main micro-mechanical parameters in the PB model, which are schematized in Fig. 6.The normal and shear stiffness of linear contact element can be respectively expressed as where A is cross-sectional area between particles, i.e.A = 2rt.The other parameters are: t = 1 (PFC 2D ), r represents the contact radius of particle-particle or particle-wall where r = min(R (1) , R (2) ) in particle-particle or r = R (1) in particle-wall, L represents the contact distance where L = R (1) + R (2) in particle-particle or L = R (1) in particle-wall.
The quantity E * represents the effective modulus of the linear contact element, and κ * represents the normal-shear stiffness ratio of the linear contact element.In addition, the linear contact element also includes the surface gap of contact g s and the frictional coefficient μ.
Similarly, the normal and shear stiffness of the parallel-bond element can be respectively written as follows where E * is the effective modulus of the parallel-bond element, and κ * is the normal-shear stiffness ratio of the parallel-bond element.In addition, σ c represents the tensile strength of the parallel-bond element, and c is the cohesion.The quantity φ is the internal friction angle, and the shear strength of the parallel-bond element reads τ c = c + σ tan φ .When the value of the stress exceeds σ c or τ c , the parallel-bond element will break and it will degenerate into the linear contact model.

Discussion on the macroscopic parameters of each single phase
Due to the lack of experimental results on each single phase, the stiffness and strength of LF and DF should be discussed firstly.Detailed mineral composition and mineral stiffness properties in these two facies were strictly confirmed by Ramos et al. 7 using SEM and CT, which are shown in Table 1.The minerals in Table 1 are quartz, calcite, dolomite, kaolinite, illite, montmorillonite, pyrite, orthoclase, and albite.The material of clay is actually the sum of kaolinite, illite, and montmorillonite.( 4)  The homogenization method based on the inclusion theory has been widely applied to predict the stiffness of rocks containing multiple mineral inclusions 41 .One useful method is the Mori-Tanaka (MT) method, which is efficiently used in the stiffness prediction of rocks with several kinds of mineral inclusions, i.e.
where C is the effective elastic stiffness tensor, C 0 represents the elastic stiffness tensor of the matrix, and N is the number of components in the shale including the matrix and (N-1) kinds of inclusions.The subscript r is the component number, where subscript 0 corresponds to the matrix.The parameter f r is the volume fraction of each component, satisfying N−1 r=0 f r = 1 .The quantity C r is the elastic stiffness tensor of the r-th component, and T r is the tensor for orientation averaging of the r-th component, i.e.
where S r is the Eshelby tensor of the r-th component, satisfying S r = P r C 0 .The expression of the Hill tensor P r of sphere inclusions can be found in the previous result 4 .
These isotropic stiffness parameters of LF and DF can be predicted by the MT method, which are demonstrated in Table 1.Overall, it can be expressed that the bulk modulus of LF is close to that of DF, while the shear modulus of LF is significantly higher than that of DF.It should be emphasized that the moduli in Table 1 are the dynamic elastic moduli 13 .The static moduli related to static mechanical properties of rock are obtained by the conversion of the dynamic elastic moduli, which are expressed as where C st represents the static elastic constant, C dyn represents the dynamic elastic constant, and α is the coefficient.The quantity α of the Young's modulus is taken as 0.39 in this study.Correspondingly, the Young's moduli of LF and DF are 33.62GPa and 26.64 GPa.
The compressive stiffness of LF and DF was determined by the uniaxial compressive strength of sandstone and soft rock obtained by He et al. 26 .The sandstone in this test contains 54% quartz, 25% clay, 13% dolomite, 2% albite, and 6% siderite.And the soft rock contains 36% quartz, 38% clay, 25% dolomite, and 1% siderite.The uniaxial compressive strength for the sandstone and the soft rock is 84.791MPa or 16.766 MPa.Compared with the composition of these two kinds of rocks, the uniaxial compressive strength ratio between LF and DF can be set to approximately 4:1.

Microscopic parameter calibration of each single phase
For the basic consensus on macroscopic and microscopic parameters of rock 9,40 5) the ratio of σ c and τ c controls the failure mode of the sample.Thus, a series of homogeneous numerical samples (25 mm × 50 mm) containing 27,222 disc particles with the porosity being 0.1 were generated to explore the relationship between macroscopic and microscopic parameters.
The microscopic parameters of the rock example in help files of PFC 2D are used as the reference parameters in this calibration process.Afterwards, four groups of numerical tests with one uniaxial tensile test and three uniaxial compressive tests were conducted in sequence for the calibration of E * , E * , к * and τ c .Each previous calibration result serves as the input parameter of next calibration test.In order to simplify the influence of various microscopic parameters, the two parameter assumptions commonly using in DEM are expressed as which means φ=0 • .
The relationships between macroscopic and microscopic parameters by fitting of LF as an example are shown in Fig. 7, and the calibration results of LF and DF are summarized in Table 2.The macroscopic and microscopic parameters take a linear correlation trend as a whole.However, it is difficult to obtain satisfactory simulation ( 6) www.nature.com/scientificreports/results directly using the calibration results for subsequent simulations, probably due to the following reasons 34,42 .The first one is the lack of the tensile Young's modulus test data, and E * may be overestimated when the tensile Young's modulus is considered equal to the compressive Young's modulus.The other one is that the uniaxial compressive strength of LF and DF still differ from that of the sandstone and soft rock, respectively, because of different content and particle size of minerals.Therefore, on the basis of calibration parameters, the revised microscopic parameters can be obtained through repeated 'trial and error' and comparison with experimental results, which are listed in Table 2.In addition, the contact properties between LF and DF are set to be the same as those of DF.

Comparison between experimental and numerical results
For this triaxial compression experiment, the final axial loading was carried out after various confining pressure loading paths 6 .Therefore, the confining pressures of 5 MPa is taken for simulation to comply with the impact of confining pressure on these experimental results.Figure 8 shows the simulation results with microscopic parameters of calibration and revision of the deviatoric stress-strain curve between experimental and numerical results.Herein, σ D is the deviatoric stress, and ε axial is the axial strain.The method for simulation of axial loading and confining pressure is by applying stress and displacement loads to the corresponding walls in PFC, which is shown in Fig. 8a.And the axial loading in simulation will stops until the current stress is below 0.7 times the peak stress.It is clearly seen that simulated stress-strain relationship in use of the revised microscopic parameters were more consistent with the experimental results than that directly using calibrated microscopic parameters, which are shown in Fig. 8b,c, respectively.Figure 8d,e shows the micro-crack distribution of revision and calibration for the numerical 45° samples.As mentioned in "Microscopic parameter calibration of each single phase" section, directly using calibrated microscopic parameters may not predict the experimental results, and there may even be significant errors.Calibrated microscopic strength parameters are too low to increase the impact of confining pressure.The simulation results of 45° and 90° samples in Fig. 8c has significant plastic deformation and a large number of cracks are generated within the DF which are almost completely discrete.
Compared with the results of 45° and 90° samples in Fig. 8b, there is a significant strength error between the numerical and experimental results in the 0° sample.Through the past experience of layered materials 24,26,39 , it can be concluded that the peak stress of 0° sample in the laminated rock with complete structures is usually not much lower than that of the 90° sample.It is clearly observed in Fig. 3 that there are initial cracks in the 0° sample, while there are no obvious initial cracks in the 45° and 90° samples.Meanwhile, it can be analyzed from Fig. 8d that tensile cracks can be clearly identified in the sample after testing, which should be used to compare with the generated cracks in the CT slices.And the differences of crack characteristics between the 0° sample and the other two samples in Fig. 2 may be caused by the initial defects in the 0° sample.
The initial main crack in the 0° sample is extracted by using the similar method in "Image processing of the shale sample slices" section, and the corresponding particles and contacts are removed within the main crack range.Figure 9 shows the 0° numerical sample with main cracks before compression, and the simulation results are compared with the experiments of the 0° sample in Fig. 8b.The initial crack significantly reduces the peak stress of the numerical sample, and the stress-strain curve is closer to the experimental results.In addition, the initial crack also changes the crack's evolution form in numerical results, and this fact can explain the generation of the secondary cracks in the 0° sample CT slices.If there are no obvious initial defects in the 0° sample, the cracking mode of the 0° sample should be close to oblique cracking along the main crack similar to the 45° and 90° samples.www.nature.com/scientificreports/Futhermore, macroscopic mechanical parameters of samples obtained by the experimental and numerical results are compiled and listed in Table 3.All macroscopic parameters of both from the experiment and simulation demonstrate the same variation trend with the layered anlge, which means the revised microscopic parameters in Table 2 can be applied to analyze the mechanical behaviors of the laminated shale sample.

Sensitivity analysis of microscopic mechanical parameters
When the discrete element model contains two different mineral combinations (LF and DF), more microscopic parameters in the laminated model are required to be analyzed than that in the homogenous model.The superposition of microscopic parameters in this dual-model increases the complexity of the relationship between macroscopic and microscopic parameters.Therefore, the five key microscopic parameters containing E * = E * , κ * = κ * , σ c , c and φ in DF and LF are set to 0.5, 0.75, 1, 1.25 and 1.5 times benchmark values of microscopic parameters in Table 2 (rev) for sensitivity analysis, where E * and к * are the microscopic stiffness parameters and σ c , c and φ are the microscopic parameters related to strength.After analyzing the stress-strain curves of a total of 150 models in 10 groups, the effects of each microscopic parameters on the macroscopic parameters (E and σ c ) of 0°, 45° and 90° sample are organized and expressed in Fig. 10.With the increasing of the microscopic parameters, the values of macroscopic E and σ c exhibit different responsive rule.
It can be found from Fig. 10a-d that microscopic stiffness parameters not only affect the macroscopic stiffness parameters E, but also have an impact on the strength parameters σ c .With the increasing E * in LF (E * LF ) shown in Fig. 10a, E of samples decreases and σ c of samples increases.Correspondingly, E and σ c of samples increase   10b.By comparing Fig. 10a,b, it can be observed that E * LF particularly affects E of the 90° sample whereas that E * DF affects E of the 0° and 45° samples more significantly.In terms of strength, it is interesting that E * significantly affects σ c of the 90° sample, while σ c of the 0° and 45° samples is also slightly affected by E * .Figure 10c,d shows the influence by the change of κ * in LF and DF, respectively.Although κ * usually affects the macroscopic Poisson's ratio of the sample, it can be observed that it also affects E and σ c of the samples.With the increase of κ * LF in Fig. 10c, E of the 90° sample decreases, σ c of the 90° sample increases, and macroscopic parameters of the other two samples have no significant change.
The effect of σ c , c and φ in LF and DF on macroscopic parameters is shown in Fig. 10e-j.From these data, it can be seen that these microscopic strength parameters have no significant impact on the macroscopic stiffness (E) of samples.In terms of strength, σ c, LF affects σ c of the 0° and 90° samples (Fig. 10e) and σ c, DF affected the σ c of 0° and 45° samples (Fig. 10f).The value of c LF affects σ c of the 0° and 90° samples and c DF affects σ c of the 0° and 45° samples, and the laws are similar, but the degree of impact is not the same.Especially, σ c of the 90° sample significantly increases with the increase of c .For the influence of φ in Fig. 10i,j, φ has an irregular effect on σ c .Although φ can increase the microscopic shear strength τ c , φ does not have a significant impact on σ c of sample as σ c and c.

Effect of the component distribution
We then discuss the influence of microstructures of laminated shale on the simulation.The CT images obtained in "Establishment of the numerical sample" section are used to establish numerical models that are closer to the actual situation than the ideal model.The microscopic parameters used for simulation are taken from Table 2.
The comparison of the CT model with the ideal model is shown in Fig. 11, and it can be seen that there is a significant difference in the stress-strain curves between the two models.Although the microscopic parameters are the same, it is surprising that the difference in microstructures has a significant impact on the numerical specimens.The anisotropic patterns of the CT model is consistent with that of the ideal model.That is, the value of σ c for the 0° sample is the biggest, and that of the 45° is the smallest.The value of E for the 90° sample is bigger than those of the other two samples, as they also have the same value.However, the anisotropic degree of the CT model is lower than that of the ideal model.For microstructures, the ideal model is the most anisotropic model, and the component arrangement in the CT model is not strictly an ideal layered structure.For ideal models of 0° and 90° samples, the through cracks during failure need to penetrate the LF components, so the strength of the ideal model is relatively high.However, the ideal model of the 45° sample is prone to generating oblique    on macroscopic mechanical properties is not as significant as that of the irregular distribution of components (CT), which means that the impact of the irregular distribution of components requires further optimization.

Effect of different confining pressures
In the experiment, the confining pressure is an important factor that affects the mechanical behaviors of shale especially in deep and ultra-deep reservoirs.As a result, the stress-strain curves of three kinds of layered angle samples under different confining pressures are obtained by changing the confining pressure of the numerical sample under revised microscopic parameters in Table 2, which are shown in Fig. 13.From the stress-strain curve of the three kinds of samples, it can be seen that the confining pressure significant increases the yielding and peak stress of the samples, but has no significant effect on the stiffness of the three samples.Figure 13 also shows the number and distribution of cracks in three sets of samples after failure of samples.It can be seen that the number of cracks in the failure sample increases with the increase of the confining pressure, and the main crack of samples changes from the oblique main crack to the 'X-shaped' main crack with the increase of the confining pressure.
In addition, there are some different patterns in the 0°, 45°, and 90° samples with the increasing confining pressure, and some unique properties only occur in certain samples.In the elastic stage of the stress-strain curve, the tangential modulus of the 0° sample exhibits a slight decrease during the loading process, which becomes apparent as the confining pressure increases.Then, after the elastic stage, the failure mode of the sample gradually changes from brittleness to ductility with the increase of the confining pressure, which is particularly evident in the 45° sample, and there are even plastic deformation features which appear in the 45° samples.Meanwhile, this plastic property is not evident in the 0° and 90° samples.After the sample reaches the peak stress and then it is destructed, the sample exhibits a certain residual strength that increases with the confining pressure.

Conclusions
In conclusion, a discrete element numerical sample is established to match the laminated shale, and the microstructures and mineral compositions of the laminated shale based on CT images provided by digital cores are fully considered in the modeling process.The mechanical behaviors of the shale sample are simulated via the PFC 2D , which are in accordance with the experimental results.The results show that, the laminated shale can be simplified as a laminated composite material containing both light facies (LF) and dark facies (DF) by observing CT images.The discrete element model with ideal layered structures can still simulate the macroscopic anisotropic mechanical behaviors of the laminated shale when the microscopic parameters of the LF and DF are fully considered.Through the analysis of the CT images and comparison of the experimental and simulated results, it was found that the initial cracks are the key to the difference between the experimental and simulated results.The Initial cracks lead to decrease the compressive strength of the 0° sample, and they can change the failure mode of the 0° sample.For the two-phase composite material when analyzing the laminated shale, the microscopic strength parameters of the parallel bonding model mainly affect these macroscopic strength parameters, while the microscopic stiffness parameters have a significant impact on both the macroscopic stiffness and strength

Figure 1 .
Figure 1.3D reconstruction of Mancos shale's digital core, where the label of sample follows the original data in Digital Rock Portal, i.e.PL and PD denote the parallel bedding and the perpendicular bedding, respectively, and pre and post denote before and after the triaxial test, respectively.

Figure 2 .
Figure 2. Sliced images of shale sample in the x-z plane, where (a) shows the process of slice, and (b-d) are the 3D digital cores and 2D slices before and after the triaxial compression test with the bedding angles being 0°, 45° and 90° samples.

Figure 3 .
Figure 3. Grayscale sliced images in the x-z plane, where (a-c) are the slices of 0°, 45° and 90° samples, respectively.

Figure 4 .
Figure 4.The morphological processes of sliced image of 45° sample.

Figure 5 .
Figure 5. Numerical model of 45° sample, where (a) is the binary CT image of the 45° sample, (b) is the numerical discrete element CT model, (c,d) are the local details of binary image and numerical model, respectively, and (e) is the ideal numerical model with similar layered thickness.

Figure 6 .
Figure 6.Schematic diagram of the linear parallel bond model.
, one has: (1) the effective modulus E * of the parallel-bond element controls the macroscopic Young's modulus under tension; (2) the effective modulus E * of the linear contact element and E * control the macroscopic Young's modulus under compression; (3) the quantities κ * and κ * control the Poisson's ratio; (4) the quantities σ c and τ c control the macroscopic compressive strength (σ c ); (

Figure 8 .
Figure 8.Comparison between calibration and revision microscopic parameters, where (a) is the loading mode of the numerical sample, (b,c) are the stress-strain curves with calibration and revision microscopic parameters, respectively, (d,e) are the number and distribution of cracks in 45° samples with calibration and revision microscopic parameters, respectively.

Figure 9 .
Figure 9. Simulation results of the 0° sample with initial main cracks, where (a) is the 0° numerical sample with initial cracks, (b) is the comparison among the experimental results, the ideal model and the crack containing model, and (c) is the failure mode of the 0° sample with initial cracks.

Figure 11 .
Figure 11.Comparison between the ideal and the CT models.

Figure 12 .
Figure 12.The influence of layered thickness on ideal models.

Figure 13 .
Figure 13.The stress-strain curve and crack number distribution of the 0°, 45° and 90° samples under different confining pressure values.

Table 1 .
, Mineral composition and stiffness parameters of LF and DF.

Table 2 .
Microscopic parameters of calibration and revision.

Table 3 .
Comparison of macroscopic parameters between experiment and simulation.